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£h ] Abstract 

The Riemann theta function is a complex-valued function of g complex variables. It appears 
in the construction of many (quasi-) periodic solutions of various equations of mathematical 
physics. In this paper, algorithms for its computation are given. First, a formula is derived 
allowing the pointwise approximation of Riemann theta functions, with arbitrary, user-specified 
^\ , precision. This formula is used to construct a uniform approximation formula, again with 

arbitrary precision. 

O 

^ ■ 1 Motivation 

O. 

An Abelian function is a 2 g- fold periodic, meromorphic function of g complex variables. Thus, 
Abelian functions are generalizations of elliptic functions to more than one variable. Just as elliptic 
functions are associated with elliptic surfaces, i.e., Riemann surfaces of genus 1, Abelian functions 
J> | are associated to Riemann surfaces of higher genus. As in the elliptic case, any Abelian function 

can be expressed as a ratio of homogeneous polynomials of an auxiliary function, the Riemann 
theta function [|llj. Many differential equations of mathematical physics have solutions that are 
written in terms of Abelian functions, and thus in terms of Riemann theta functions (see ||, [?J 
and references therein). Thus, to compute these solutions, one needs to compute either Abelian 
functions or theta functions. Whittaker and Watson |p0| , §21.1] state "When it is desired to obtain 
definite numerical results in problems involving Elliptic functions, the calculations are most simply 
performed with the aid of certain auxiliary functions known as Theta- functions." Whittaker and 
Watson are referring to Jacobi's theta functions, but the same can be said for Abelian functions, 
using the Riemann theta function. This paper addresses the problem of computing values of the 
Riemann theta function and its derivatives. 
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2 Definition 



The Riemann Theta function is denned by 

o(z\n)= e 2 ^ n - an+n - z ), (1) 

where z G C 9 , f2 G C 9 * 9 , such that fi is symmetric (fi T = f2) and the imaginary part of f2, Im(fi), 
is strictly positive definite. Such an fi is called a Riemann matrix. Also, nz = Yli=i n i z ii the scalar 
product of the integer vector n with z; and n ■ fi • n = j=i QijiHrij. The positive definiteness of 
Im(f2) guarantees the convergence of (|l]), for all values of z. Then the series (|l|) converges absolutely 
in both z and fi, and uniformly on compact sets. Thus, it defines a holomorphic function of both 
z and fi. The function defined by (|l]) is also known as a multidimensional theta function, or 
as a theta function of several variables. There are as many different conventions for writing the 
Riemann theta function as there are names for it. These different conventions differ from @) by at 
worst a complex scaling transformation on the arguments. An extensive overview of the wealth of 



properties of 0(z|fi) is found in [14, [15], 16] 



Remarks 

• The Riemann theta function was devised by Riemann as a generalization of Jacobi's theta 



functions of one variable 12], for solving the Jacobi inversion problem on general compact 



connected Riemann surfaces [17]. For these purposes Riemann considered only theta functions 
associated with Riemann surfaces. Riemann theta functions in their full generality, as defined 
in (|l]), were considered first by Wirtinger [21], whose convention for the arguments is adopted 



here. 



Often, so-called Riemann theta functions with characteristics are considered [14]. Such theta 
functions with characteristics are up to an exponential factor Riemann theta functions (|l]) 
evaluated at a shifted argument. Thus the computation of the Riemann theta function (|l]) 
also allows the computation of theta functions with arbitrary characteristics. 

In many applications, the Riemann theta function ([!]) originates from a specific Riemann 
surface, i.e., the Riemann matrix fi is the normalized periodmatrix of the Riemann surface. 
Obtaining this periodmatrix is a nontrivial problem, which is now also reduced to a black-box 
program. This issue was addressed in ||. 

Some of the algorithms given in this paper have been used already by some of the authors (For 
instance, tori with constant mean curvature, and Willmore tori with umbilic lines were con- 
structed in [|l0| using the results of ||, |j ; multiphase solutions of the Kadomtsev-Petviashvili 
equation were constructed using Schottky uniformization in [||, |5)) but they were not easily 
available for use by others. Now, these algorithms are implemented as black-box programs 
in Maple and Java. The maple implementation is included in the Maple distribution as 



of Maple 8. It is also available from [http : //www . math . f su. edu/^hoei j/RiemannTheta/ 



The Java implementation is available from www-sfb288.math.tu-berlin.de/~jem. These 
implementations are discussed in the appendices. 
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3 Rewriting the Riemann theta function 

The Fourier series representation @) is the starting point for the computation of Riemann theta 
functions. Separating both z and ft in their real and imaginary parts, z = x + iy, ft = X + iY, 
and using Y = Y T , we obtain 



6(z\Sl) = V e Mh n - an + n ' z ) 

nezn 

_ e 2Tri(±n-X-n+n-x+i(±n-Y-n+n-y)) 

ueZb 

_ e ^i(^n-X-n+n-x) e -2Tr(^n-Y-n+n-y) 

e 2 m (\n-X-n+n-x) ^{lin+Y^yyY ■(n+Y- 1 y)+-\y-Y- 1 -y) 

y ^ e 2ni(±n-X-n+n-x) e -K(n+Y- 1 y)-Y-(n+Y- 1 y) ^ 



At this point, all exponential growth has been isolated: the factor multiplying the sum grows 
double-exponentially as the components of z leave the real line, due to the fact that Y~ l is strictly 
positive definite: this follows since Y is strictly positive definite, which also guarantees the existence 
of Y^ 1 . All terms remaining in the sum are either oscillating (the first factor of every term), or a 
damped exponential (the second factor of every term) . 

Since for most applications, the exponential growth term cancels out, it will be disregarded in 
almost all that follows: any statements of approximation, pointwise or uniform, of the 
Riemann theta function pertain to the infinite sum in (|^), without considering the 
exponential growth. 

Let [V] be the vector with integer component closest to V, and let [[V]] = V — [V]. Continuing 
the rewriting of the theta function, 

6{z\ft) = e^Y- 1 *^ e ^(^.X.n + n. a; ) e -K^[^-V] + [[^-V]])^{r l+ [Y- 1 y] + [[Y- 1 y]]) 
= e ,y Y-"y £ e 2. i (i(n-[r- 1 2/ ]).X.(n-[r- 1 2/ ]) + (n-[r- 1 2/ ]). a; ) x 

xe - 7r (n + [[Y- y ]]).F.(n + [[F-y]])_ (3) 

The last step is achieved by shifting the summation index n. From this formulation, it is clear that 
the size of each one of the terms in the infinite sum is controlled by its second exponential factor. 

Similarly, formulae are worked out for the derivatives of theta functions. Denote the iV-th order 
directional derivative of the g- variable Riemann theta function along the ^-dimensional vectors , 
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Then 



= (27Ti) N e^- Y ' 1 -y J2 (k (N) • {n - [Y^y])) . . . (fc (JV) • (n - [Y~ l y 



x 



xe - 7r (n + [[^ y ]]).F.(n + [[F- y ]]) ) (g) 

using similar steps as before. The exponential growth is still factored out, but the remaining infinite 
sum now contains terms that grow algebraically as the argument of the Riemann theta function 
leaves the real line. The order of this growth is equal to the order of the derivative. However, 
the size of individual terms within the infinite sum is again determined by the last factor which is 
exponentially decaying. 

Remark. Riemann theta functions of genus greater than one have been computed before, for 
instance in ||, where genus three was considered. There four distinct representations of these Rie- 
mann theta functions were used, based on whether the different phases behaved as trigonometric or 
hyperbolic functions (two limits of elliptic functions), depending on the parameters in the Riemann 
matrix. These distinct representations were required to have a manageable number of contributing 
terms to the series. Such a variety of representations is not needed here: all limit cases are incor- 
porated in the form (|3|) . This is obtained by the separation of the real and imaginary parts of both 
the argument z and the Riemann matrix Q,. 



4 Pointwise Approximation 

The formulae @ and @ are the basis for the approximation of theta functions and their derivatives. 
All approximations are based on determining which terms of the infinite sum are dominant. 
For the Riemann theta function, 



x e 



Sr 

x e -„(n + [[Y-iy)]).Y{n + [[Y-y])) ^ 



where 

S R ={ne Z> (n + [[Y-'y]]) ■ Y • (n+ [[Y^y]}) < R 2 } . (7) 

Below, we show this limit exists, hence proving that the Riemann theta function is well defined. 
This proof is different from that found in many references (see [14], for example) in that it also 
allows the estimation of the error made in the approximation of the truncation of the series by only 
considering a finite value of R. 
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Since Y is strictly positive definite, it has a Cholesky decomposition, Y = T T T. Let A be 
the lattice of all vectors v(n) of the form v(n) = y^vrT (n + [f^" 1 ?/]])! f° r n G ^ s '• Then is 
rewritten as 



S fl = |«(n) € A 



l«(n)|| < 



}■ 



(8) 



Thus, the approximation requires finding all lattice points in A that are also inside the (^-dimensional 

sphere | \v(n)\ \ = R. 

Let 

e R (z\n) = e -y^"-y^ e ^<K n -[ y " 1 ^])- x < n -[ Y " 1 ^]) + ( n -[ r "^])- :E ) e -ii«wii 2 , (9 ) 



Sr 



then 



e(R) 



\e( z \n) -e R (z\n)\e-*y- Y K y 

^ e ^(K n -[^"^])- x -( n -[ r "^) + ( ri -[ r "^)- a; ) e -ii^) 



(10) 



A\S R 



^(i(n-[r- 1 y]).X.(n-[r- 1 y]) + (n-[r- 1 ?/ ]). a; ) || W(n) 



e -||«(n) 

a\s r 



(11) 



Thus e(i?) is the absolute error of the oscillatory part of the Riemann theta function, due to 
the truncation to a finite number of terms. In practice, this oscillatory part is of order 1, so e(R) 
is also a measure of the relative error. 

In order to estimate this last sum, the following theorems and lemmas are used. 

Theorem 1 (Mean- value Theorem for subharmonic functions) Let Q, be a domain in M. 9 . 
Let u be a twice continuously differ entiable function on £1, continuous on the boundary of d£l oftt, 
which is subharmonic on f2. Then for any ball B R (y) in Q of radius R and center y 



u(y) < 



1 



u(x)dx. 



(12) 



B^VolBf (0) J B 9 r(v) 
A proof of this theorem is found in || . 

Lemma 1 For every integer p > and any dimension g > the function u{y) = e~"2/" 2 ||?/|| p with 
y £ M. 9 is subharmonic on M. 9 \ B R (0), with R = \ \j g + 2p + \f~g 2 ~ + 8p. 

Proof. Since u(y) depends only on r = \\y\\, this calculation is done in ^-dimensional spherical 
coordinates. 



. , . d 2 u a — 1 du 



2u{y) 2r l - (g + 2p) + 



p 2 + (g - 2)p 
2r 2 



This last factor is positive if r > 9 + 2p + \J g 2 + 8p, from which the result follows. 
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Lemma 2 Let A be a g- dimensional affine lattice in R 9 , i.e., A = {Xn + x\n G Z 9 }, with X G 
G7(n,R), x G R 9 , and p G Z ; positive. Then 

£ lli/ll^ ll! " p <f(-) ! 'r(^,(«-p/2f), (is) 



where T(z, al) = t z 1 e z dz, the incomplete Gamma function j^j, R > ^ y g + 2p + y^ 2 + 8p + 
p/2, and /? = min {| \x — y\ \ \x, y G A, x ^ y}. 

Proof. By the previous lemmas, 

V I l-^l If e -M2/ll 2 < V - I \\x\\ v e-\\ x W 2 dx 

yeA,\\y\\>R yeA,\\y\\>R yH/ ' iy ' JB P /2 [ y> 

1 (2 



V y [ \\x\\Ve-\\ x W 2 dx (14) 



< , (1Y [ ll^lPe-H^I 2 ^ 



VolSf(O) \pj jR a \B a R _ p/2 (o) 

Voisf(o) Vp7 

g fiy r e - H(9+P) /2-i dt 

2 VP/ J(R-p/2)2 
9 ( 2 V + P 



, . . , -,(R-p/2)' 

2 \pj V 2 

Here VolSr^OVVolT^O) = 5 /r was used. ■ 

Remark. The inequality following (|l4|) is quite crude. Specifically, it adds all the space in 
between the balls around the lattice points. Taking this in to account, another estimate is as 
follows: 

2 1 ( 2\ 9 f ,, „„ _||-w||2 



V Ilt/ll p e-ll^l < - - V / \\x\\ p e-^dx 

yeA,\\y\\>R 1V ; VF/ yeA,\\y\\>R B p/2^y> 

2 Vc(A,p/2) / \\x\\P e -^ 2 dx. 

J Jk9\ B ° r _ p/2 {o) 



VoLBf(O) \p 



Here c(A, p/2) is the "fill factor" of the lattice A with balls of radius p/2: it is the fraction of 
the lattice volume that is filled if a ball of radius p/2 is placed at every lattice point. This esti- 
mate is justified, since the function being integrated only depends on the radial variable. Thus, 
within thin radial shells the function is almost constant, independent of the proximity to a lat- 
tice point. The crudeness of the inequality following (|ll|) is illustrated in Fig. [I], where the 
fill factor c(A, p/2) is shown for the scenario of minimal |A| = p 9 , i.e., a square lattice. Then 
c(A,p/2) = Tr 9 / 2 /{2 9 - 1 gT(g/2)), which decays fast as g —* oo. Figure |l| shows that for a square 
lattice, not using the fill factor results in the error being overestimated by only one digit for g = 6. 
Using the fill factor for generic lattices, where A > p 9 , usually leads to a more significant improve- 
ment. 
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Thus c(A,p/2) = VolB»(0)/|A| 



Figure 1: Fill factor c(A, p/2) for the case of a square lattice. 
2vrf/ 2 (p/2)9 



-, where |A| is the determinant of the lattice. It 



p/2^/i-i 5 r (5/2 ) | A | 

is the volume of a cell of A spanned by a set of generating vectors. This gives 



l/eA,||!/||>fl 



| y ||p e -nvii" £ voigf" 1 ^) W {pW 
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VoiBf(o) Vp7 sr( 5 /2) |A| y*_ p/2 

^ 1 1 r- 



-r 2 r 9-l+ P(ir , 



5r( 5 /2) |A| 2 . 

^9/2 



|A|rfo/2) 

Using (|l3|) with p = 0, (|TTj) becomes 



r(^,(i?-p/2) 2 



e(R)< J; e -ll«WH a < | (^r(|(E-p/2f), 



(15) 



(16) 



where p is the length of the shortest lattice vector in A: p = min{||x||,x S A}, and R > (y / 2^ + p)/2 
(to satisfy the subharmonicity condition of Lemma [l]). As claimed, this proves that the Riemann 
theta function is well-defined, since lim^oo T(z,R) = Q. 

Thus, in order to approximate the oscillatory part of a Riemann theta function of g complex 
variables with a pre-determined error e, one solves the equation 



, (fl - p/2) 2 ) 



(17) 

2g + p)/2 suffices. These 



for R, with R > + p)/2, real. If no solution exists, then R = 

results combine to give 

Theorem 2 (Pointwise Approximation) The Riemann theta function is approximated by 



Sr 



■ x ) e -\\v{n)\\ 2 



(18) 
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with absolute error e on the sum. Here Sr = {v(n) G A| ||u(n)|| < R}, A = 
|y / 7rT(n + \n £ Z 9 }. The shortest distance between any two points of A is denoted by 

p. Then the radius R is determined as the greater of (y/2g + p)/2 and the real positive solution of 
e = g2 9 ~ 1 T{g/2,(R-p/2) 2 )/p 9 . 

This theorem gives a pointwise approximation to the Riemann theta function: for every z at 
which the Riemann theta function is approximated, Sr is different, although R is unchanged as 
long as g and e remain the same. The error used in Theorem |2| is referred to as the 100% Error 
(100%E). 

Another good estimate for R is obtained from (|l5|). Then R is the greater of (\[2g + p)/2 and 
the real positive solution of e = tt^ 2 T (g/2, (R - p/2) 2 ) /(\A\T(g/2)). This error is referred to as 
the Fill Factor Error (FFE). 

Figure ^ compares the size R of the ellipsoid, required to obtain a prescribed error e, as a result 
of using the 100%E or the FFE, for the cases g = 2 and g = 16. In these figures, the worst-case 
scenario |A| = p 9 was assumed. As expected, the difference between the 100%E and the FFE is 
small for small genus, but becomes significant for larger genus. 




(a) g = 2 (b) g = 16 

Figure 2: R, assuming the 100% Error (solid line) and Fill Factor Error (dashed line) as functions 
of — logio(e), which is the number of accurate digits. 

Two examples are shown in Table [l|, for various values of e. Both of these use the Riemann 
matrix Q, g with Qjj = i, j = 1, . . . , g, and Qjk = —1/2, j ^ k. The first example (left side of the 
table) corresponds to g = 2, the second (right side) to g = 6. For both examples, the Riemann theta 
function is evaluated at z = 0, thus in this case the oscillatory part of the Riemann theta function 
equals the Riemann theta function. The value of the Riemann theta function is given in the first 
line. The first column of an example gives e, used for the pointwise approximation, with either the 
100%E, or the FFE. In the table N(100%E) and N(FFE) denote the number of terms that are used 
to compute the oscillatory part of the Riemann theta function, i.e., the number of elements of 5r, 
using the 100%E or the FFE, respectively. AE denotes Actual Error, i.e., the difference between 
the actual value of the oscillatory part of the Riemann theta function (computed using 30 digits of 
accuracy, and e = 10~ 30 ) and the computed value, using the two different error formulae. 

Using (^), we see that the first example, with g = 2, computes n2 =-oo cos(nin27r)e _7r ^ rai+n2 \ 
Since for this example g = 2, there should be little difference between the computation using 
the 100%E or the FFE. This is confirmed: although the value of R is slightly different for both 
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computations, the number of elements of Sr is identical, resulting in both computations having 
the same accuracy. 

For the second example, with g = 6, there is a difference between the two computations, 
although it is not as outspoken as for g = 16 in Fig. |2[ It is clear that the Actual Error using the 
FFE computation is closer to e, than using the 100%E. Even in this case, the computation does 
several orders of magnitude better than prescribed. This, of course, is due to the inequalities that 
were used to obtain these error estimates. 





9 = 


2 








g = 6 






6>(0,0|ft 2 )e-^- r 2 V 


= 1.1654010572 


0(o,o,o,o,o,o|n 6 ) e -^- r 6 y 


= 8.3721839831 


e 


iV(100%E) 
= TV(FFE) 


AE(100%FE) 
=AE(FFE) 


e 


iV(100%E) 


AE(100%FE) 


iV(FFE) 


AE(FFE) 


1E-1 


5 




7.5E-3 


1E-1 


485 


4.3E-5 


233 


9.2E-4 


1E-2 


9 




1.5E-5 


1E-2 


797 


3.8E-6 


485 


4.3E-5 


1E-3 


13 




1.2E-6 


1E-3 


1341 


2.6E-7 


797 


3.8E-6 


1E-4 


21 




5.1E-11 


1E-4 


2301 


1.3E-8 


1341 


2.6E-7 


1E-5 


21 




5.1E-11 


1E-5 


3321 


4.2E-10 


2301 


1.3E-8 


1E-6 


21 




5.1E-11 


1E-6 


4197 


3.8E-11 


3321 


4.2E-10 


1E-7 


21 




5.1E-11 


1E-7 


5757 


2.3E-12 


4197 


3.8E-11 


1E-8 


25 




1.9E-12 


1E-8 


8157 


9.2E-14 


5757 


2.3E-12 


1E-9 


29 




1.8E-13 


1E-9 


10237 


3.5E-15 


8157 


9.2E-14 


1E-10 


37 




1.5E-17 


1E-10 


12277 


2.7E-16 


10237 


3.5E-15 



Table 1: Two examples of using the pointwise approximation to compute the oscillatory part of a 
Riemann theta function. 



Remarks 

• The method for approximating the oscillatory part of a Riemann theta function requires the 
determination of the elements of Sr. This is the set of all elements of A that lie inside the g- 
dimensional sphere determined by ||w(n)|| = R. It is easier to consider the equivalent problem 
of determining the points with integer coordinates that lie inside the ^-dimensional ellipsoid 
determined by 7r(n — c) • Y • (n — c) = R 2 , where c is the center of the ellipsoid. This is 
done recursively, by remarking that every {g — l)-dimensional plane section of a g-dimensional 
ellipsoid is a (g — l)-dimensional sphere. A finite number (due to the discrete nature of the 
lattice) of such (g — l)-dimensional sections are taken. All of these are {g — l)-dimensional 
ellipsoids, on which this section process is repeated, until one arrives at a one-dimensional 
ellipsoid, i.e., a line piece. At this level it is easy to determine which integer points are in this 
piece. This method takes full advantage of the triangular form of the Cholesky decomposition 
of Y = T T T. This is done as follows: all integer points satisfying 

\\T(n-c)\\<R/V^ (19) 
are sought for. Since T is upper triangular, this implies \T gg {n g — c g )\ < R/t/n, or 

R R 
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This gives a set of allowed values of n g . For each one of these, write n = (n, n g ) T , c = (c, c g 
With 

t i 

T 

u ± gg 

19| ) becomes ||T(n — c) + — c 9 )|| 2 + T 2 g (n g — c g ) 2 < R 2 /tt, which is rewritten as 



\T(h -c + T t(n g - cg))\\ < JR 2 /tt - T$Jn g - c g Y 



This is a similar problem as (|I9|), but in dimension g — 1 instead of g, and with a different R, 
c and T. The procedure is now repeated. 

Both error estimates require the knowledge of p, the shortest lattice vector of A. The problem 
of finding the shortest lattice vector of a given lattice is one of the main problems in the study 
of lattices. It is usually addressed using lattice reduction. Lattice reduction attempts to find 
a standard form for the matrix of generating vectors tj, j = 1, . . . ,g of the lattice, in which 
certain intrinsic (i.e., independent of the representation of the generating vectors; such as p, 
for instance) properties of the lattice are easier to examine. The problem of computing p is 
known to be NP hard in the parameters g and ln(maxj(| \tj\\)) ||i~9|| . However, an approximate 
algorithm due to Lenstra, Lenstra and Lovasz (the LLL algorithm, |{13[|) is known which 
is polynomially hard in g and ln(maxj(||ij||)). This approximate algorithm is exact in low 
dimensions (guaranteed in 1,2 and 3), but only gives approximate answers in high dimensions. 
In practice, the algorithm has been found very satisfactory. A rigorous error estimation is 
possible using the results of the LLL algorithm, as it provides both a lower and an upper 
bound for p: If p is the value of p found by the LLL algorithm, then p/2^ 9 ^ 1 ^ 2 < p < p. 
Since the LLL algorithm for relatively low genera (g < 10) usually finds the correct value of 
p, this will not be pursued here. 



5 Uniform Approximation 

It is now easy to extend the results of the previous section to obtain results for the uniform 
approximation of the oscillatory part of Riemann theta functions. In Theorem 2, the sum extends 
over all terms whose summation index vector is inside an ellipsoid. The size of this ellipsoid (R) 
is determined by the allowed error of the pointwise approximation. The shape of the ellipsoid 
depends on CI, but not on z. The center of the ellipsoid is — [[y _1 y]]. So, if different arguments z 
are considered, only the center of the ellipsoid changes and not its shape or size. But, this center 
can only wander over a (/-dimensional cube of volume 1, centered around the origin. This leads to 
the following 

Theorem 3 (Uniform Approximation) The Riemann theta function is approximated by 

^|fi)«e^ r "^e 2 ^KM r ^M^ (20) 

Ur 

with absolute error e on the sum. The approximation of the sum is uniform in z. Here v(n) = 
^T{n+[[Y- l y]]), 

U R = {n£ Z 9 \ir(n - c) Y ■ (n - c) < R 2 ,\cj\ < 1/2, j = l,...,g} . (21) 
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Let A = {^JirTfo + [\Y 1 y]}) \n E Z 9 }. The shortest distance between any two points of A is 
denoted by p. Then the radius R is determined as the greater of (y/2g + p)/2 and the real positive 
solution ofe = g 2V~ l T (g/2, (R - p/2) 2 ) / p9 . 

Thus, to obtain a uniform in z approximation with absolute error e for the oscillatory part of 
the Riemann theta function parametrized by fi, it suffices to add all terms whose summation index 
is inside the envelope of all ellipsoids of size R (determined by e) whose shape is determined by ft, 
and whose center is at any point inside the unit cube centered at the origin. This object is also 
obtained as the convex hull of all ellipsoids of size R placed at all corners of the unit cube centered 
at the origin. This situation is illustrated in Fig. ||. The ellipsoid on the left illustrates the scenario 
of obtaining the summation indices for a pointwise approximation. The deformed ellipsoid on the 
right illustrates the case of a uniform approximation. The deformation is a consequence of the 
wandering of the center of the ellipsoid over the unit cube. 



o o 



o o 



it o 



it II 



it o 



it o 



it o 



it o 



(a) 



(b) 



Figure 3: Summation indices (a) for a pointwise approximation, inside an ellipsoid centered at 
[[Y --1 ?/]], and (b) for a uniform approximation, inside an ellipsoid whose center moves around the 
unit cube. 

It is clear that using a uniform approximation results in more terms of the sum being used. 
In many applications (graphics, creation of value tables, etc.) many values of the same Riemann 
theta function (i.e., with constant ft) are computed. In such cases it is beneficial to use the same 
representation for the computation. 



Remarks 

• For the purposes of obtaining a uniform approximation formula, it is essential that the expo- 
nential growth was factored out in (|3|). 

• The above theorem uses the 100 % Error of Theorem ^[ Just as in the case of a pointwise 
approximation, it is possible to consider a uniform approximation using a Fill Factor Error. 
This results in similar modifications to Theorem || as for Theorem |2[ 

• Determining the set Ur is not trivial. In practice, it may be more convenient to work with a 
single ellipsoid centered at the origin, which is too large. 

As an example, consider the Riemann matrix 
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f 1.690983006 + 0.9510565162 i 1.5 + 0.363271264 i \ 

~\ 1.5 + 0.363271264 % 1.309016994 + 0.9510565162 i )' ^ ' 

This is the Riemann matrix associated with the genus 2 Riemann surface obtained by compactifying 
and desingularizing the algebraic curve fi — X 7 + 2A 3 /i = 0. It was computed using the algorithms 
outlined in ||. Using the above theorem with e = 0.001, 23 terms have a contribution to the 
oscillatory part of the Riemann theta function. They have summation indices 

U R (W0%E) = {(-2, 2), (2, -2), (0, 2), (0, -2), (2, -1), (-2, 1), (2, 0), (-2, 0), (0, -1), 
(0, 1), (-1, -1), (1, 1), (1, 0), (-1, 0), (-2, -1), (2, 1), (-1, 1), (1, -1), 
(0, 0), (-1, -2), (1, 2), (-1, 2), (1, -2)} . (23) 

Thus a sum of 23 terms suffices to approximate the oscillatory part of the Riemann theta function 



associated with (22) with an absolute error of e = 0.001. This approximation is valid for all 
z = (x± + iyi,X2 + W2)- These 23 terms should be compared to the 12-17 terms required for a 
pointwise approximation. The pointwise approximation uses sets of summation indices which are 



subsets of (23), but these subsets differ depending on the value of z = {x\ + iyi,X2 + £2/2)- 

The uniform approximation can be used to graph various slices of the oscillatory part of the 
Riemann theta function. This is a complex function of four real variables, of which many different 
slices are possible. Twelve graphs are shown in Fig. ||, using the uniform approximation. Note 
that the oscillatory part of the Riemann theta function is periodic with period 1 in the real part of 
all components of z. This periodicity is inherited from the Riemann theta function. 

6 Derivatives of Riemann theta functions 

In this section, approximation results for derivatives of Riemann theta functions are obtained. The 
results are not as clean as for the Riemann theta function itself: even after removing the exponential 
growth, algebraic growth terms remain in the sum. The order of the algebraic growth equals the 
order of the derivative. There is only growth on the imaginary axes of the arguments z. Therefore, 
approximations of the sum with absolute error are bound to be valid only pointwise in z. It is 
possible to obtain uniform approximations in z which are valid in a bounded area of C 9 . We give 
the details for directional derivatives of first and second order. Results for higher-order directional 
derivatives follow in a similar fashion. 

6.1 First-order derivatives 

With N = 1 and = k, f| gives 

e-rV'Y^-V D(k)9(z\n) 
= 2mk-Y, (»- [Y-'V]) e 2 -<K"-[^])^<"-[ y "^]) + K[^])) e -II^WII 2 



Xe -\Mn)\\\ (24) 
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(a 4 ) (b 4 ) (c 4 ) 



Figure 4: Various views of the oscillatory part of the Riemann theta function parametrized by the 
Riemann matrix given in (^), with FFE=0.001. All plots are oscillatory parts of 6{x + iy, 0|f2) 
(index 1), 0(0, x + iy\fl) (index 2), 6(x,y\Q) (index 3), 8(ix,iy\il) (index 4). Shown are the real 
part (a), the imaginary part (b) and the absolute value (c). 
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where Sr and v(n) are defined as in (||). As before, let 
D R (k)e(z\Q,) = e*V :Y ~ K V x 

x2aifc£ (n-[r- 1 y ])e 2 ^<K n 4^ 1 ^)- X < n 4^ 1 ^) + K[ r ^])) e -ll^)ll 2 ,(25) 



Sr 



then 
e(R) 



\D(k)8(z\n) - D R (k)6(z\n)\e 



(26) 



< 



2^ fc .^(n-[r"^])^<K-[^])^(--r^) + Kr^)) e -ll^) 
27r||fc|| \\ n ~ [^"^]|| e" l|U(n)l12 



A\S f 



2TT\\k\\ 2 ^ T_1 ^H " [[^]] " [^^] 



A\S fl 

2vr||fc|| ^2 

A\S R 



\v(n) 



)=T- l v{n) - Y^y 

ITT 



\v(n) 



—T- l v(n) 

TT 



< 2ir\\k\\ V 

A\S R 

< 2V7r||fe|| (IT" 1 !! ^ ||«( 



n e 



A\S, 



e HI"(n)ll 2 + 27 r||fe||||r- 1 y|| J] e -ll^)H 2 

A\S R 

- M ^)H 2 +27r||fe||||y- 1 j,|| J] e-II^WH 2 , 

A\S R 



(27) 



where the Cauchy-Schwarz and triangle inequalities were used. The second term in ( |27D is estimated 
as before. For the sum in the first term, Lemma [l] is used with p = 1. However, this only applied 
in the region where | |i)(n)| |e~" t '^ n ^' is subharmonic. According to Theorem [l], this holds when 



(28) 



R > \ v 9 + 2 + vV + 8 - Thus 

'2V 



e (i?) < v^3l|fe|| (llr- 1 1| r (£±i, (a - p/2) 2 ) + ^ ||y-i y || r (|, (n - p/2) 2 )) , 

provided i? > \ (^Jg + 2 + y / g 2 + 8 + ^ . 

Theorem 4 (Pointwise Approximation of the First Directional Derivative) The direc- 
tional derivative of the Riemann theta function along k is approximated by 



D{k)e{z\n) = e^y- Y l y x 

x2vrzfc £ (n- [Y^y]) e ^<| (r*- [V"-^] ) -Jt- (r*- [^V] [l-- 1 ?/] )) e _| 1^(^)1 1- ? (29) 

Sr 

with absolute error e on the scalar product of liri k with the sum. Here Sr = 
{v(n) e A\ \\v(n)\\ < R}, A = {^T(n + [[Y^y]}) \n G Z 9 }. The shortest distance be- 
tween any two points of A is denoted by p. Then the radius R is determined 
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as the greater of ly g + 2 + y / g 2 + 8 + pj /2 and the real positive solution of e = 

y/* g \\k\\ (l) 9 (\\t-'\\ r (R - p/2f) + ||y-i y || r (1 , (R - p/2f)). 

This theorem clearly gives a pointwise approximation: not only is the location of the ellipsoid 
Sr dependent on the evaluation point z; more importantly its size R depends on it. This is the 
reason why this result cannot be extended to a uniform in z approximation theorem, as in the case 
of Theorem ||. By limiting the domain of y, it is possible to obtain a uniform in z approximation 
theorem, valid over this domain: 

Theorem 5 (Uniform Approximation of the First Directional Derivative) The direc- 
tional derivative of the Riemann theta function along k is approximated by 



D(k)6(z\n) = e^y- Y l y x 

x2nik-Y / (n-[F- 1 y ])e^<K n 4^ 1 ^)- X < n -[ F ^]) + ( n -[^ 1 ^)) e -ll^)ll 2 ,(30) 

Ur 

with absolute error e on the scalar product of2irik with the sum, where the approximation is uniform 
in z, z = x + iy, for y G -B£(0). Here v(n) = v / 7rT(n + [[V -1 ?/]]), and 

U R = {ne Z 9 \ir(n-c) ■ Y • (n - c) < R 2 ,\cj\ < 1/2, j = l,...,g} . (31) 

Let A = {■ v /7rT(n + [[l^ -1 ?/]]) |nsZ 9 }. The shortest distance between any two 
points of A is denoted by p. Then the radius R is determined as the greater of 

^g + 2 + y / g 2 + § + p^J /2 and the real positive solution of e = y/irg\\k\\ ||T _1 || 
r (2±i, (R - p/2f) + \\t-'\\ r (1 , (R - p/2f)). 

Proof. This follows easily from the Pointwise Approximation result, by using ||l^ _1 y || < 1| \\y\\ 

and WY^W < WT^P. ■ 



This theorem is not as useful as Theorem ^ due to the restriction on the size of When 
it is used, it results in a situation where many terms are included which are only relevant for the 
evaluation of values of the directional derivative near \\y\\ = L. However, for graphing or other 
purposes where many evaluations of D(k)0(z\ft) are required, centered around an area for z where 
y = 0, such a uniform approximation is usually beneficial compared to the pointwise approximation. 
It was used in Fig.JBa to illustrate the linear growth of the derivative of the Riemann theta function 
parametrized by (|22), after removal of the exponential growth. 

6.2 Second-order derivatives 

With N = 2, (§) gives 

e^y- Y ^-y D(k^\k^)e(z\n) 

= (27Ttf • (»- [ Y ^y])) ( fc(2) • (»- [ Y ^y 
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20000 
1 0000 


-1 0000 
-20000 




(a) 



(b) 



Figure 5: Two directional derivatives of the Riemann theta function 9(0,ix\Q) with Q as in (p2|). 
As before, the exponential growth-factor has been factored out. In (a) the directional derivative 
of order one with k = (1,0) is shown. Fig. (b) illustrates the quadratic growth of the directional 
derivative of order two, with k^ 1 ' = = (1,0). 

xe2m (i(n-[F- 1 y]).X.(n-[l-- 1 y]) + (n-[F- 1 I/ ])) eH |^ (n) || 2 
= (2™)Mim {n-lY^y])) (k^ . („_ [Y^y])) x 



Sr 



x e 



2m (i(n-[F- 1 y]).X.(n-[r- 1 2/ ]) + (n-[r- 1 ?/ ])) eH | V(n) 



(32) 



where Sr and v(n) are defined as in (||). Using similar steps as before, and using Lemmas [l] and 
^ with p = 2 results in the following theorems: 

Theorem 6 (Pointwise Approximation of the Second Directional Derivative) The 

directional derivative of second order along k^ and k^ of the Riemann theta function is 
approximated by 



D(kW,kW)d(z\n) = e*y- Y ~ 1 -v x 

Sr 

xe 2 7r i(i(n-[F- 1 y]).X.(n-[Y- 1 2/ ])+(n-[y- 1 y])) e _|| V(n) 



(33) 



with absolute error e on the product of (27ri) 2 with the sum. Here Sr = {v(n) G A| ||u(n)|| < R}, 
A = { yfixT(n + [[Y ~ 1 y]]) |n£Z 9 }. The shortest distance between any two points of A is de- 
noted by p. Then the radius R is determined as the greater of f^J g + 4 + \J g 1 + 16 + p \ /2 

and the real positive solution of e = 27r^||fc (1) || ||fc (2) || UY (pT^T (R - p/2) 2 ) + 



2^||T- 1 || Hr-^H r (2±i, (R- p/2) 2 + vr l^yf T (§, (R - p/2) 2 
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Theorem 7 (Uniform Approximation of the Second Directional Derivative) The directional 
derivative of second order along k^ and k^ of the Riemann theta function is approximated by 

D(k^,k^)9(z\n) = e ^y- Y ' 1 -y x 

Ur 

xe2 .<i(n-[r- 1 y]).X.(n-[r- 1 y ]) + (n-[r- 1 2/ ])) e _ || , (n)||2 {M) 

with absolute error e on the scalar product of (2iri) 2 with the sum, where the approximation is 
uniform in z, z = x + iy, for y E S£(0). Here v(n) = yfirT(n + [[l^ _1 y]]), and 

U R = {n E Z 9 \ir(n-c) ■ Y • (n - c) < i? 2 ,|cj| < 1/2, j = 1,... , 5 } . (35) 

Lei A = {■ v /7rT(n+ [["K _1 y]]) |n E Z 9 }. T/ie shortest distance between any two points of A is 

denoted by p. Then the radius R is determined as the greater of fyj g + 4 + \fg^ + 16 + /2 

and toe real positive solution of e = 2vrc/||fc (1) || ||fc (2) || fjQ || T_1 || 2 ( r ( £ ^. ~ Z 5 / 2 ) 2 ) + 2 V^ 
WT-^LT^^R- p/2) 2 ) + TT^HT^frd^E-p^) 2 ) 



The same remarks as for the first-order directional derivative are valid here as well. The last 
theorem has limited use due to its restriction on the size of ||y||. It was used in Fig. ||b to illustrate 
the quadratic growth of the derivative of the Riemann theta function parametrized by (^), after 
removal of the exponential growth. 



7 Siegel transformations 

Consider the Riemann theta function parametrized by the Riemann matrix 

-1 / 111.207 96.616 
~ 2iri V 96.616 83.943 

This is the example from Appendix C of p|, adapted to the definition of the theta function used 
here. It was used in || to illustrate the need for a fundamental region of Riemann matrices. 
The problem arising is that the ellipsoid (ellipse, in this case) determining the summation indices 
in the approximation of the oscillatory part of the Riemann theta function is very eccentric. The 
eigenvalues of Y are 31.0587 and 0.000324, resulting in an eccentricity of the ellipse of 1—0.54 10~ 10 . 
Thus very few of the summation indices closest to (0, 0) play a part in the evaluation of the Riemann 
theta function. This was the problem addressed in ||. Since our algorithm incorporates a way to 
determine which summation indices lie inside the ellipsoid determined by Y, this is not troublesome 
here. What is troublesome is that the ellipsoid contains many integer points: already for e = 0.001, 
the evaluation of #(0,0|f2) requires the inclusion of 109 terms. This is largely due to the fact that 
the ellipsoid lies along a rational direction in the (ni,rt2) plane, as illustrated in Fig. |(| 

In such cases, the transformation properties of the Riemann theta function can be used to great 
advantage. Just as the elliptic functions and the Jacobian theta functions, the Riemann theta 
function has a modular transformation property. The modular transformation is a transformation 
on the Riemann matrix. The transformed Riemann matrix defines a Riemann theta function, just 
as the original Riemann matrix does. These two Riemann theta functions are related: up to an 
affine transformation of the argument z and an overall scaling factor, they are identical. For details, 
see Q. 
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(a) 



(b) 



Figure 6: The summation indices for the evaluation of the oscillatory part of 0(0, 0|f2) and 6(0, 0\ft) 
with e = 0.001. In Fig. (a), with ft, there are 109 summation indices, all inside a very eccentric 
ellipse lying along the line ri2 = —n\. Figure (b) corresponds to ft, which was obtained from ft 
by way of a modular transformation (|37|). Now only one summation index lies inside the ellipse: 
(ni,n 2 ) = (0,0). 

Remarks 

• Geometrically, the modular transformation on the Riemann matrix is a transformation on 
the period lattice of the Riemann theta function. Thus the modular transformation property 
relates two Riemann theta functions with different period lattices flZf , 

• For Riemann matrices originating from Riemann surfaces, the modular transformation amounts 
to choosing a different canonical intersection basis for the homology of the Riemann surface 
|7j. Then it is not a surprise that the solutions of differential equations written in terms of 
Abelian functions do not depend on this choice. One could think of this as a discrete symme- 
try of the solutions of the differential equation. Thus for such applications it is unimportant 
to transform the Riemann theta function, as long as one consistently does all calculations 
with a preferred homology basis. 

Let 

be a symplectic matrix with integer elements, T € SP(2g, Z), i.e., 

( a b \ ( g I 9 \(a T c T \ = (0 g I g \ 

\ C (1 J \ ~Ig Og J \ b T (F J \~I g Og )> 

where a, b, c, and d are g x g matrices with integer elements, I g and g are the g x g identity and 
nul matrix, respectively. SP(2g, Z) is called the modular group. An element of the modular group 
transforms the Riemann matrix ft according to 

ft -> ft = (aft + b)(cft + dy 1 . (37) 
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The modular group is generated by the three generators Ti,^ and T3 [14]: 
fa 



g la 



l g b 



% lg 

Og "I 

■9 % 



9 nT ), n = afta \ 6(Az\AQA T ) = 9(z\ft), (38) 
, fl = fl + b, 9(z\Cl + b) = #(z + diag(b)/2|ft), (39) 



r 3 =(7 9 ), n = -n- 1 , e(n- l z\ - n- 1 ) = y/det(-m) e wiz - n - z e(z\n), (40) 



where diag(b) denotes the diagonal part of the matrix b, which has even diagonal elements, and 
b T = b. In the last equation, the branch of the square root is used which is positive when its 
argument is positive. 

Since any modular transformation is a composition of these three types, the transformation 
formulas on the Riemann theta function on the right can be used to calculate the effect of the 
resulting modular transformation on the Riemann theta function. 

From our point of view, the goal is to find a transformation on the Riemann matrix (i.e., a 
composition of generating modular transformations) to minimize the eccentricity of the ellipsoid 
determined by the transformed matrix. The algorithm described in this section is due to Siegel 
]l~8| . Siegel's goal was to construct a fundamental region for Riemann matrices, analogous to the 
elliptic case. The algorithm iteratively finds a new Riemann matrix with improved (i.e., smaller) 
eccentricity for the ellipsoid it determines. However, the algorithm is not optimal. An optimal 
algorithm does not appear to be known. The description below is an algorithmic description of 
Siegel's [18], aimed at implementation. This implementation is included in Maple 8. 

Theorem 8 (Siegel Reduction) Every Riemann matrix fi can, by means of a modular trans- 
formation $3%), be reduced to a Riemann matrix ft = X + iY,Y = T T T, with X (Y ) the real 
(imaginary) part of CI, and T upper triangular, such that 

1. \X jk \<l/2, j,k = l...g, 

2. the length of the shortest lattice vector p of the lattice generated by the columns of T is bound 
from below by <J V3/2, and 

3. max{\Nj\ : ||TiV|| < R,R > 0, fixed, N G Z 9 } has an upper bound which only depends on g 
and R. Thus this upper bound is independent ofT. 



The second condition eliminates ellipsoids with high eccentricity. The theorem in effect guarantees 
an upper bound for the number of summation indices required for the evaluation of the oscillatory 
part of a genus g Riemann theta function with a prescribed error e (and thus R). This upper bound 
does not depend on the reduced Riemann matrix ft. 

Proof: 

• Proof of statement 1: Using (|39|), 

n-> «- [Re(«)]. (41) 
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Proof of statement 2: We can assume that T is lattice reduced, using the LLL algorithm 
1 13]. Then p = T\\ = a/Yh, since T is upper triangular. Siegel [18] shows that the determi- 
nants of the imaginary parts of two Riemann matrices connected by a modular transformation 
( |37| ) are related by 

where Y, Y are the imaginary parts of f2 and ft, respectively. Using the modular transfor- 
mation with 

T \ , / -1 T \ / 1 T \ , / T 



with the (g — l)-dimensional zero vector, ( |42| ) becomes 

|det(y)| 



|det(r)| 



l^lll 2 



If [On I < 1, this transformation is applied. This transformation preceeded by the trans- 
formation giving © is repeated until |On| > 1. Then |fin| 2 = X 2 X + Fu > 1. Thus 



p = vYii > y y 1 — > Y v / 3/2, since X 2 X < 1/4. It remains to be shown that this 
iteration terminates. This is demonstrated by Siegel (l8|] . 

• Proof of statement 3: If T is lattice reduced using the LLL algorithm, then (see fl~3|| ) 

* Tjj > 0, 1<3<9, 

* \T jk \ < \Tjj\/2, l<j<k<g, 

■ 7 :",-i • 7 / i,-: T Jr >' •/• <■>■ 
It follows from these last two properties that Tjj > ^-Tj-x^-i, for j G (l,g], and thus 
2ij > Mf\ 3 1 p. Then max{|A^| : ||TJV|| < R,R > 0, fixed, N G Z»} < max{|iy,-| : 
|(TJV)j| <R,R> 0, fixed, j e[l,g],N € 'L 9 } < J2/r. ■ 

Note that the preceding proof indeed shows that the new ellipsoid defined by [|TJV|| = R has an 
eccentricity bound away from 1. 

This algorithm was used on the Riemann matrix given in (|36|), with great success: after applying 
the algorithm, the evaluation of the oscillatory part of the Riemann theta function parametrized 
by the transformed Riemann matrix at z = (0, 0) with absolute error e = 0.001 requires only one 
summation index: (ni,n2) = (0,0). The modular transformation used is 

f 8 7 \ 
7 6 
6-700 
\ -7 8 / 

resulting in 

A _ / 7.94597 -3.94937 
~ V -3.94937 14.4545 

giving an ellipsoid with eccentricity 0.927921, with a major-axes ratio ~ 3/8. 
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Appendix A: The Maple implementation 

The maple code can be viewed by typing the following commands in Maple 8 or later: 

interface (verboseproc=2) ; 

op(RiemannTheta) ; 

op ( ' RiemannTheta/ doit ' ) ; 

op( 'RiemannTheta/boundingellipsoid' ) ; 

op ( ' RiemannTheta/ f indvectors ' ) ; 

op('RiemannTheta/make_proc') ; 

op( 'RiemannTheta/f initesum' ) ; 



It can also be downloaded from [http : //www . math . f su . edu/^hoei j /RiemannTheta/ 



The procedure RiemannTheta uses a variety of arguments. The first two of these are required. 
All others are optional and can be omitted. 

1. A g x g Riemann matrix ft. 

2. A g-dimensional vector z. 

3. A number e indicating the desired accuracy of the result. 

4. A (possibly empty) list of vectors. Each vector is used as in (Q) for the calculation of directional 
derivatives. 

5. Other optional arguments: The algorithm computes the exponential growth factor a and the 
oscillatory part b separately. The user can specify if a and b should both be given in the 
output, or if the output should consist of only exp(a)6. 

RiemannTheta distinguishes two cases: 

a. z evaluates to an element of <C 9 . 

b. z does not evaluate to an element of <C 9 , because it contains one or more variables that have 
no assigned values. 

If no directional derivatives are given, then RiemannTheta computes either a list of two complex 
numbers [a, b] or one complex number e a b such that 0(0,, z) ~ e a b. Here b is bounded for z G C 9 , 
and has error less than e. If directional derivatives are given, then b grows polynomially. For case 
b, it is not guaranteed that 6's error is bounded by e unless an a-priori bound for z is known. 

If RiemannTheta is called, then the set S of points in the "ellipsoid" described in Figure || 
is calculated. See Figure ||(a) for case a, and Figure |(b) for case b. Following the computa- 
tion, RiemannTheta returns an expression that contains 'RiemannTheta/f initesum' (args). Here 
'RiemannTheta/f initesum' is a procedure and args is a list of arguments. These arguments con- 
tain all information necessary to compute the value of RiemannTheta with the specified accuracy. 
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This information includes z, fi and S, but reorganized into a form more convenient for summation 
(e.g., the real and imaginary parts of the input are separated, as are the integer and non-integer 
part before summation over the set S is begun). However, the summation is not done unless z £ (D 9 '. 
Only if z contains no unassigned variables can z be evaluated to an element of C 9 , and only in 
this case will 'RiemannTheta/f initesum' (args) return a number. This number then becomes 
the output of RiemannTheta. If z does contain unassigned variables, then the summation will be 
delayed, and 'RiemannTheta/f initesum' (args) will remain unevaluated. Thus, using variables 
in z instead of complex numbers causes RiemannTheta to return an answer that contains an un- 
evaluated procedure. This is useful if one wants to calculate more than one value of RiemannTheta 
for the same fi. If only a single value of RiemannTheta is wanted, one gives RiemannTheta a vector 
z £ (D 9 . If one wants to compute RiemannTheta for many z's (e.g., for plotting purposes) then 
one includes variables in z. The procedure which is the output of RiemannTheta is then used for 
evaluation. This way the set of points in Figure |3](b) is computed only once, which is more efficient 
than computing the points in Figure ||(a) for each value of z. 

Below is a verbatim example of interactive use of Maple and RiemannTheta. 



22 



> r3 := sqrt(-3)/3: 

> M := evalf( Matrix (2 ,2 , [ [l+2*r3, -l-r3] , [-l-r3 , l+2*r3] ] ) ); 

[ 1. + 1.154700539 I -1. - 0.5773502693 I] 

M : = [ ] 
[-1. - 0.5773502693 I 1. + 1.154700539 I ] 

> eps := 0.001: 

> z : = [1-1, 1+1] ; # Here z contains no variables. 

z := [1 - I, 1 + I] 

> R := RiemannTheta(z , M, [] , eps); 

-8 

R := -21.76547256 - 0.2323298326 10 I 

> L := RiemannTheta(z , M, [] , eps, output =1 i st ) ; 

-10 

L := [3.627598727, -0.5785248137 - 0.6175311504 10 I] 

> a := L[l] : 

> b := L[2] : 

> R := exp(a)*b; 

-8 

R := -21.76547256 - 0.2323298326 10 I 

> z := [X, 1+1] ; # Now z contains 1 variable. 

z := [X, 1 + I] 

> R := RiemannTheta(z , M, [] , eps, output =1 i st ) ; 

2 

R := [3.627598726 Im(X) + 3.627598724 Im(X) + 3.627598726, 
RiemannTheta/f initesum( . . . ) ] 

> eval(R, X=l-I) ; 

-10 

[3.627598728, -0.57852736323 - 0.61753198638 10 I] 

> eval(R, X=1-2*I) ; 

-9 

[10.88279618, 0.62464131574 - 0.56009279234 10 I] 

> eval(R, X=1-3*I) ; 

-10 

[25.39319109, 0.44006321314 - 0.92901122197 10 I] 



Appendix B: The Java implementation 

The Java implementation is an adaptation of an earlier implementation in C which has been in use 
since 1994. The Java version uses the definition ([l]) of the Riemann theta function. The C version 
employed a different definition: 

9 (z\n) = y e \n B. n +z n = q ( * *\ 
y 1 ' ^ \ 2m 2m J 

where the real part of B is negative definite. 

The java implementation of the Riemann theta function is realized in the package riemann . theta, 
which includes the public classes LatticePointsInEllipsoid, ModularGroup, SiegelReduction, 
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bias 



mfc 



11 



nemann 



numericalMethods 



11 



riemann.theta 



+ Theta 

+ ThetaWithChar 

+ TransformationPropertySupport 

+ ModularPropertySupport 

+ SiegelReduction 

+ ModularGroup 

+ LatticePointsInEllipsoid 



Figure 7: Dependencies of the package riemann.theta. 



ModularPropertySupport, TransformationPropertySupport, Theta and ThetaWithChar. Unlike 
other programming languages, Java does not incorporate a complex type, thus the package requires 
an implementation of such, realized in mfc . number . Complex. The implementation also requires the 
packages bias (basic linear algebra system), and numericalMethods. This last package is used only 
in internal computations. Figure ^ illustrates these dependencies, showing only the classes within 
riemann.theta. 

The class Theta implements the uniform approximation theorems |3[ |5| and |7[ The default 
error used by the program is the FFE, with default error tolerance set to 10 . Siegel's reduction 
algorithm as described in section ^ is always used. 

Example: Consider the slices of a genus 2 Riemann theta function shown in Figure ||. The 
complete Java program for computing the data in one of these slices is: 

import riemann. theta. Theta; 

import mfc. number. Complex; 

import blas.ComplexVector; 

import blas.ComplexMatrix; 

public class ExampleSlice { 

public static void main( String [] argv ) { 

// create the period matrix 

ComplexMatrix B = new ComplexMatrix( 2 ); 

// set the period matrix 

B.set( 0,0, 1.690983006, 0.9510565162 ); 

B.set( 0,1, 1.500000000, 0.3632712640 ); 

B.set( 1,0, 1.500000000, 0.3632712640 ); 

B.set( 1,1, 1.309016994, 0.9510565162 ); 
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} 



// use the different normalization 

B.setTimes( new Complex ( 0, 2 * Math. PI ) ); 

// create a complex argument vector 

ComplexVector V = new ComplexVector ( 2 ); 

// create a theta function instance with period matrix 

Theta theta = new Theta( B ) ; 

// create storage for the slice 

Complex [] [] slice = new Complex [101] [101] ; 

// loop over grid 

for( int i=0; i<101; i++ ) 

for( int j=0; j<101; j++ ) { 
// set grid vector 

V.set( 1, 2 * i *Math.PI / 100, 2 * j * Math. PI / 100 ) ; 
// evaluate theta function at grid vector 
slice [i] [j] = theta. theta( V ); 

} 

} 



To compile this example, the user needs the Java archives riemann. jar, bias, jar, mfc. jar and 
numericalMethods . jar in their classpath. These archives, as well as their source code and docu- 
mentation is available from www-sfb288.math.tu-berlin.de/~jem. 

The package riemann. theta also includes a class ThetaWithChar which incorporates Riemann 
theta functions with characteristics. 
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